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Abstract. The entrainment transition of coupled random frequency oscillators 
presents a long-standing problem in nonlinear physics. The onset of entrainment in 
populations of large but finite size exhibits strong sensitivity to fluctuations in the 
oscillator density at the synchronizing frequency. This is the source for the unusual 
values assumed by the correlation size exponent v' . Locally coupled oscillators on a d- 
dimensional lattice exhibit two types of frequency entrainment: symmetry-breaking at 
d > 4, and aggregation of compact synchronized domains in three and four dimensions. 
Various critical properties of the transition are well captured by finite-size scaling 
relations with simple yet unconventional exponent values. 
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1. Introduction 

Mathematical modeling of coherent oscillation in a large population of autonomous units 
has a long and interesting history, as discussed in several monographs [TJ [2] and recent 
reviews [H HI |5]. One particular paradigm is the Kuramoto model of coupled random 
frequency oscillators [2J. In terms of the phase variables <pj,j = 1, . . . , N, the dynamical 
equations take the form, 

^^^-KAsm^-e). (1) 

Here K is the coupling strength among the oscillators, and A(> 0) and 6 are the 
amplitude and phase of the mean instantaneous phase factor, 

N 

Ae ie = N^^e^', (2) 
i=i 

which is also known as the order parameter for entrainment. The intrinsic frequencies 
Uj are drawn from a distribution g(u). In this paper we shall focus on the case where 
g(uj) is a gaussian function centered at w max = with unit variance. 

In the classic work by Kuramoto[2j, existence of an entrained state is established 
by considering solutions to ([1]) at a constant A. Under this assumption, Eq. ([1]) is 
easily integrated. After an initial transient, oscillators with \ojj\ < KA each reaches 
a fixed angle and together form the entrained group, while those with \ujj\ > KA go 

into periodic motion at modified frequencies Uj = — (KA) 2 . In the limit N — > <x>, 
only the entrained population contributes to the sum (j2]). The resulting self-consistent 
condition 

/KA 
dug{uo)yJl - (u/KA) 2 (3) 
-KA 

has a nontrival solution A > when K > K c = 2/irg(0). Thus entrainment requires the 
product Kg(0) to exceed a threshold value. On the supercritical side, A(K) ~ (K—K C )P 
with the usual mean- field exponent /3 = 1/2. 

While the above predictions agree well with numerical solutions of the globally 
coupled Kuramoto model, there are several unresolved issues regarding the stability of 
the mean-field solution, as discussed elegantly in Ref. [3]. From a physicist's point of 
view, one would like to understand how temporal fluctuations of A(t), which are present 
when the sum ((2]) contains only a finite number of terms, affect the behavior predicted by 
the mean-field treatment. When the oscillator frequencies Uj are drawn independently 
from the distribution g(u), sample-to-sample fluctuations in the time averaged A need 
to be considered as well. These effects culminate at K = K c where the size of the 
entrained group grows sub-linearly with the system size N. 

Daido developed a perturbation theory[6] to calculate the susceptibility defined by, 

x = N(& - A 2 ) (4) 

Here the overline bar denotes time average. His calculation yielded a power-law 
divergence of x a ^ K c , but the exponent 7 differs on the two sides of the transition: 
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7 = 7_ = 1 on the subcritical side and 7 = 7+ = 1/4 on the supercritical side. Daido 
has further suggested a finite-size scaling form based on his own numerical studies, 



Here k = (K — K c )/K c is the distance to the transition, v'_ = 2 and u' + = 1/2. These 
exponent values unfortunately do not agree with later numerical work[7]. 

Equations (0Q) have been adapted to describe synchronization phenomena in finite 
dimensions with local oscillator coupling [8], 



Here Aj denotes the set of the nearest neighbor sites of j on a d- dimensional hypercubic 
lattice. Arguments have been presented to show that, as the linear system size L — > 00, 
global entrainment is not possible when d < 2[8j. Analytical and numerical studies of 
the model in three to six dimensions by Hong, Chate, Park and Tang (HCPT)[9] have 
established the following. For d > 4, the entrainment transition is of the symmetry- 
breaking type with exponents indicative of a mean-field behavior. At d — 3 and 4, on the 
other hand, global entrainment is reached through aggregation of locally synchronized 
domains as the coupling strength K increases beyond a critical value K c . Although 
the aggregation process is continuous and can be associated with a length scale £ that 
diverges at K c , the fraction of entrained oscillators may undergo a discontinuous jump 
at K c in the limit L — > 00. 

Below we shall give an overview of our current knowledge on the finite-size scaling 
properties in the all-to-all coupled Kuramoto model and discuss some fine differences 
between this case and the Kuramoto model on a random graph. We shall argue that 
the latter provides the correct mean-field description of the entrainment transition at 
d > 4. Key features of the entrainment process in three and four dimensions are also 
discussed. However, no attempt is made to review the vast amount of related work in 
the literature. Fortunately, this role is partially fulfilled by several recent articles cited 
above. 

2. Finite-size scaling in the globally coupled model 

2.1. Anomalous behavior due to quenched frequency fluctuations 

In typical numerical studies of ([I]), the frequencies Uj are drawn independently from the 
distribution g(u). HCPT have shown that, in this case, v' = 5/2 on either side of the 
transition. Since the argument is quite straightforward and relevant for the discussion 
of synchronization on complex networks and in finite dimensions, we reproduce it below. 

Still assuming A in (p]) to be time-independent and restricting the sum in ([2j) to 
the entrained oscillators in a finite population, Eq. ([3]) is replaced by 




(5) 




(6) 




j,\ojj\<KA 



(7) 
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Since the actual values of the ufs vary from sample to sample, the sum \1/ has a 

1 /2 

"quenched" fluctuation 5^/ = \E r — (\&) oc N s /N, where N s is the number of oscillators in 
the frequency interval (— KA, KA). Here (•) denotes average over different realizations 
of the random frequencies. Close to the transition, an expansion of Eq. ([?]) at small A 
yields, 

A = (K/K C )A - c(KA) 3 + 5V, (8) 

where c = —7ig"(0)/lQ. Mathematically, the variance of 5^ can be computed from 
(ED for independently drawn frequencies: ((5$) 2 ) = 4g(0)KA/3N + 0(A 2 /N). Thus 
solution to OH]) takes the scaling form, 

A(K,N) = N-^fikN 2 / 5 ). (9) 

The scaling function f(x) is sample-dependent and satisfies the equation, 

xf - cK\f + ^/^f^tf 1 ' 2 = 0. (10) 

Here e = 5^> / ((5\1 / ) 2 ) 1//2 is a gaussian random variable with zero mean and unit variance. 

Equation ([9]) suggests v' = 5/2 as seen in numerical experiments. Its origin lies 
in the fluctuations in the oscillator density at the peak frequency cj max of g(co). This 
effect produces a sample-dependent shift 5K C ~ jV~ 2 / 5 of the entrainment threshold. At 
the nominal threshold K = K c , A ~ 5K@ ~ N^ 1 ^ 5 . In other words, the number of 
entrained oscillators right at K = K c scales as N 4 ^ 5 . (Note that, when the term e in 
( [TO]) is negative, there is only the trivial solution / = at x = 0. For these samples, the 
entrained population size is much smaller than A^ 4//5 .) Numerical simulations of (TTDT) 
produced values of (A 2 ) in quantitative agreement with that of the Kuramoto model 
for K > K c $\. 

2.2. Dynamic fluctuations in the subcritical region 

For K < K c , oscillators essentially run with their intrinsic frequencies. If the phase 
factors in the sum (EJ) were all statistically independent from each other, one would 

2 

have A = 1/N and hence x — !• The perturbative calculation by Daido|6j yielded the 
following expression for K close to K c : 

x = — - — + D. (11) 
X K c -K K ' 

Here A = 4/K c for the gaussian g(u), and D is a numerical constant which can also be 
evaluated, though the calculation is rather tedious. Equation (TTTj) was rederived recent 
by Hildebrand et al. using a different approach |10j. 

Numerical results by Hong et a/.|ll] are in good agreement with (llip when the 
system size > N c ~ \K — K c \~ 5 / 2 . At smaller values of N, a sample-dependent 
behavior is seen, so the result depends on the averaging procedure. For example, 
the data for \ calculated from (jlj) for a given sample and then averaged over many 
samples obey the scaling (x) — |A;| _1 \l/(/cA r2/ ' 5 ). This implies that, in the critical region 
\k\ < N~ 2 / 5 , dynamic fluctuations of A have an amplitude 5Ad yn ~ N~ 3 / 10 , which is 
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weaker than the sample to sample fluctuation A ~ iV -1 / 5 . Due to this difference, the 
hyperscaling relation v' = 2/3 + 7 is violated. At present, there is no analytic theory 
that explains these numerical findings. 

2.3. Finite-size scaling in the "pure" case 

The oscillator density fluctuations on the frequency axis responsible for the anomalous 
v' = 5/2 can be eliminated when the frequency Uj of the j'th oscillator is chosen 
according to the formula 



The intrinsic frequencies of the oscillators in this case are nearly uniformly spaced locally 
on the axis. In analogy with the discussion of disordered systems, we call this situation 
the "pure" case. 

Numerical simulations of the Kuramoto model with frequencies given by (j!2p 
produced exactly the same A (if) function in the large size limit as the "disordered" 
case. On the other hand, measurement of fluctuations of A(t) yielded j± = 1/4 and 
u' ± = 5/4 which are very different from the disordered case [11]. It thus appears that the 
subcritical behavior ( II lj) requires extra assumptions which are not fulfilled in the pure 
case. 

2.4- Entrainment transition on complex networks 

The entrainment transition of the Kuramoto model has been studied on complex 
networks by several groups j5]. Here the dynamical equations are given by (jSJ), with 
Aj denoting the set of nodes connected to node j on the network. In the case of scale- 
free networks where the degree k of nodes satisfies a power-law distribution P(k) ~ k~ a , 
both numerical and analytical studies show that the exponent (3 = 1/2 for a > 5, but 
increases sharply as = l/(a — 3) for 3 < a < 5 [12]. The size of the critical region where 
finite-size effects become significant, as described by the exponent z/, also broadens in 
the latter case. These two effects combine to give a broadened transition when the local 
connectivity of the network assumes a broad distribution. 

Interestingly, for a > 5, and on the random graph with a fixed degree, the finite- 
size scaling form ([5]) is found to hold across the transition, with the exponents given 
by 7± = 3/2 and v' ± = 5/2 [13]. This is in contrast to the behavior of the Kuramoto 
model with all-to-all coupling. To better illustrate the difference, let us examine the 
distribution of the mean phase velocity Vj = dipj/dt of oscillators in the neighborhood of 
the transition. Figure 1(a) shows a comparison of the cumulative distributions obtained 
from the all-to-all coupling case and on a random graph at k = 6. Both systems are 
in the critical region of their respective entrainment transitions. Entrained oscillators 
have the same phase velocity t> syn ~ as indicated by the vertical jumps of the curves. 
Distinct behavior is seen right above and below t> syn in the two cases. For all-to-all 
coupling, the curve is nearly flat, suggesting a vanishing density of oscillators in the 




(12) 
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Figure 1. (Color online) (a) The cumulative distribution of the phase velocity at 
criticality for all-to-all (black, K — and random-graph (red, K = 0.38, degree 

= 6) couplings. Here N — 3200 in both cases, (b) The phase velocity distribution in 
a given sample in three dimensions at different coupling strength K. Here N = 128 3 
and K c ~ 0.7. Dashed line indicates the frequency distribution at K = 0. 

neighborhood of v syn . Thus the entrained population is well-separated from the running 
oscillators. On the random graph, the vertical jump joins smoothly with the detrained 
regions. The rise at small \v\ can be fitted to a \v | 1//2 -law. Consequently, the oscillator 
density near t> syn diverges as \v \~ l l 2 . The boundary between the entrained and running 
populations is much less well-defined in this case. 

The difference in the density of oscillators near the entrainment frequency can be 
appreciated from the following consideration. For all-to-all coupling, the fluctuation 
effect on a given oscillator is due to temporal variations in the system-averaged order 
parameter A which vanishes in the large size limit. In addition, from the scaling 
Adyn ~ jV~ 3 / 10 and Ao ~ N~ 1 ^ mentioned above, A can be regarded as a constant 
in the first approximation even in the critical region. The mean phase velocity of 
oscillators then essentially follows the behavior predicted by Kuramoto's mean-field 
analysis, which is what the black curve in Fig. 1(a) suggests. In contrast, on the 
random graph, each oscillator is coupled to a small number of other oscillators. Hence 
the fluctuation effect has a finite strength and does not decrease significantly in the 
infinite size limit. Oscillators with the "right environment" first become entrained, but 
due to the presence of fluctuations, they undergo diffusive phase motion in a weak 
ordering field A, punctuated with phase slips. This diffusive background is present even 
on the entrained side. 
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3. The Kuramoto model in finite dimensions 

The "quenched" random frequencies in (J6j) introduce a static phase deformation <fr^ 
whose spatial structure has been argued to be responsible for the different types of 
entrainment behavior in finite dimensions [9]. Salient features of this static phase field 
are captured by a linear analysis carried out by Sakaguchi et al. [8] and by Hong et al. [H] 
The resulting dynamic equations, which reduces to the "quenched Edwards- Wilkinson 
equation" [15] in the continuum limit, exhibit a diffusive relaxation towards a "fully 
entrained" state with vanishing phase velocity for all oscillators. As the linear system 
size L — > oo, fluctuations of remain bounded when the system dimension d > 4, but 
diverge when d < 4. In addition, the gradient of <ft^ diverges when d < 2. 

HCPT proposed a plausible scenario for the entrainment transition based on the 
linear analysis. The linear approximation is justified at sufficiently large K and above the 
lower critical dimension d\ = 2, where the run-away oscillators are absent. (We assume 
here that the intrinsic frequencies are bounded.) As the coupling strength K decreases, 
more and more oscillators break away from entrainment. The question is whether the 
remaining oscillators can stay phase-locked as their fraction in the population diminishes. 
For d > 4, the entrained state has bounded phase fluctuations and hence breaks the 
phase symmetry (f> — > (f> + c. In this case, A as defined by Eq. (T5j) can still be used 
as a global order parameter and acts as a spontaneous ordering field on the oscillators 
as in the globally coupled case. On the other hand, for d < 4, the phase symmetry is 
not broken. Therefore phase-locking need to be maintained via a percolative network 
of entrained oscillators in space. The latter requires a finite density of these oscillators 
given the random assignment of the oscillator frequencies. The entrainment transition, 
when viewed in terms of the fraction of entrained oscillators in the system, is expected 
to be discontinuous. 

The above picture is largely confirmed by extensive simulations carried out by 
HCPT in three to six dimensions which also yielded additional details of the entrainment 
transition [9]. In five and six dimensions, the distribution of the mean phase velocity 
shows essentially the same behavior as the Kuramoto model on a random graph as 
depicted in Fig. 1(a). On the other hand, a very different distribution of phase velocities 
is obtained in three and four dimensions. Figure 1(b) presents an example from a 
iV = 128 3 system in three dimensions. Significant narrowing of the distribution takes 
place at K values well below the threshold K c ~ 0.7 at this size. The wings of the 
distribution can be fitted to a power-law P(v) ~ M~ 2 - The data is consistent with 
the picture that, below K c , entrainment takes place locally through the formation of 
entrained clusters. The size of the clusters grow as K increases, reaching the system 
size at K c . The increase of the entrainment threshold with the system size, as observed 
in simulations, can be taken as further confirmation of this picture[9]. 

HCPT also examined various singular behavior of the system at the entrainment 
transition with the help of finite-size scaling analysis. For d > 4, A as defined by 
Eq. 02]) is a suitable order parameter. For d = 3 and 4, static phase fluctuations in 
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a sufficiently large system render A vanish even on the entrained side. Instead, the 
Edwards- Anderson order parameter 

A EA = hm l|Wto(*)-fc(*>)] (13) 
J 

was found to be appropriate. Numerical results for a range of system sizes are in good 
agreement with the finite-size scaling 

(A 2 ), (A| A ) = L~ w/u ^(kL 1/u ). (14) 

The exponent values extracted from the simulation data can be expressed as, 

d<A: (3 = 0, u = 2/(d-2); 

d > 4 : 0= 1/2, v = 5/(2d), v' = dv = 5/2. (15) 

Note that in three and four dimensions, /3 = is consistent with a jump in the fraction 
of entrained oscillators at the transition. 



4. Conclusions and outlook 



Numerical investigations of the Kuramoto model with all-to-all coupling, on scale-free 
networks and random graphs, and on finite dimensional lattices, revealed a wealth of 
behavior at the onset of global entrainment. The usual finite-size scaling approach has 
been shown to be successful in capturing the key critical properties of these systems, 
including fluctuation effects and the size of the entrained population at the transition. 
Development of phenomenological arguments based on various approximations enhanced 
our understanding of the entrainment process. However, a full analytic theory of 
fluctuations is still lacking, even in the case of the Kuramoto model with all-to-all 
coupling. 

The original Kuramoto model (CQ) has a small parameter A in the neighborhood 
of the transition which can be employed for systematic expansions. However, such 
a treatment need to be handled with care when considering oscillators with ufs 
comparable to A. These oscillators are in fact the ones that become entrained first 
and contribute significantly to A. The good news is that, as Fig. 1(a) shows, temporal 
fluctuations of A are weak when compared to A itself at sufficiently large N. Therefore 
one may try to adapt Daido's scheme j6] on the supercritical side to a finite system at 
K = K c to compute the size of the entrained population, from which the exponent v' 
can be obtained. This work is in progress. 

For the Kuramoto model on random graph and on finite-dimensional lattices with 
d > 4, the above perturbative scheme can not be applied directly. However, since the 
entrainment transition is of the symmetry breaking type, a coarse-grained description 
is expected to be appropriate. For example, one may consider equations of the complex 
Ginzburg-Landau type with suitable thermal and quenched disorder terms that mimic 
the effects of the random frequencies. Standard field theoretic techniques can then be 
applied to perform the relevant calculations. 
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Existing analytic methods may not be able to treat the entrainment transition in 
three and four dimensions. To describe in detail the breakdown of entrained clusters 
as the coupling strength K weakens, one needs to know the dynamic behavior of 
defects (grain boundaries and dislocations) that mediate phase slips between neighboring 
entrained domains. Establishment of relevant scaling relations will help one to sort out 
various types of complex spatiotemporal correlations, from which a better understanding 
of the entrainment process can be expected. 
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